package scu.maqiang.cte;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.fes.FES2T32;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class CylinderResultT3_SC2 {

	public static void main(String[] args) {
	    
	    int lineNum = 100;
	    double dis = 1.0 / lineNum;

	    double kxx = 1.0, density = 1.0, specficHeat = 1.0;
	    double E = 2.0, Nu = 0.3, alpha = 0.05;
	    double beta = E / (1 - 2 * Nu)  * alpha;

	    double theta = 0.5, endTime = 10, dt = endTime / 10000;
	    double time = 0.0, val;
	    double omega = 0.25;

	    int[] label = new int[1];
	    int[] heatLabel = {3};
	    int[] dispLabel1 = {2, 4};
	    int[] dispLabel2 = {1};
	    
	    int verbosity = 1; 
	    String dir = "CylinderResultT3\\";
	    String meshFile = dir + "MyCircle142DT.msh";
	    String heatOutputFile, XYStressYFile;
	    String temp00File;

		Mesh2T3 mesh2D = new Mesh2T3();
		mesh2D.readMesh(MeshFormat.msh, meshFile);
		int outerLabel = 50, leftLabel = 70, bottomLabel = 60;
		
		ScalarFunc outerBoundary = (xy, llabel, param) -> (xy[0] > 0.00001 && xy[1] > 0.00001) ? 1.0 : -1.0;
		ScalarFunc bottomBoundary = (xy, llabel, param) -> xy[1] < 0.00001 ? 1.0 : -1.0;
		ScalarFunc leftBoundary = (xy, llabel, param) -> xy[0] < 0.00001 ? 1.0 : -1.0;
		mesh2D.setBoundaryLabel(outerBoundary, null, outerLabel);
		mesh2D.setBoundaryLabel(bottomBoundary, null, bottomLabel);
		mesh2D.setBoundaryLabel(leftBoundary, null, leftLabel);
		
	    //�����¶ȳ���Ԫ�ռ�, ��װ��Ӧ�������
	    FES2T31 fs = new FES2T31(mesh2D);
	    SRMatrix LMatrix = new SRMatrix(fs.GetNdof());
	    SRMatrix RMatrix = new SRMatrix(fs.GetNdof());
	    ScalarFunc kxxFunc = (xy, llabel, param) -> kxx;
	    ScalarFunc rhoFunc = (xy, llabel, param) -> density * specficHeat;
	    fs.assembleThetaScheme(new ScalarFunc[]{kxxFunc, rhoFunc}, 
	    		                   new double[] {dt, theta}, 
	    		                   new BVPType[] {BVPType.COMMON, BVPType.CONSIST_MASS},
	    		                   LMatrix, RMatrix);
	    double[] x = new double[fs.GetNdof()];
	    double[] RHS = new double[fs.GetNdof()];

	    //���õ��Գ�����Ԫ�ռ�, ��װ��Ӧ�������
	    FES2T32 fs2 = new FES2T32(mesh2D);
	    SRMatrix K = new SRMatrix(fs2.GetNdof());
	    SRMatrix M = new SRMatrix(fs2.GetNdof());
	    double[] d = new double[fs2.GetNdof()];
	    double[] d0 = new double[fs2.GetNdof()];
	    double[] vel = new double[fs2.GetNdof()];
	    double[] vel0 = new double[fs2.GetNdof()];
	    double[] acc = new double[fs2.GetNdof()];
	    double[] acc0 = new double[fs2.GetNdof()];
	    double[] heatLoad = new double[fs2.GetNdof()];
	    double[] totalLoad = new double[fs2.GetNdof()];
	    double[] u = new double[fs.GetNdof()];
	    double[] v = new double[fs.GetNdof()];
	    
	    fs2.assembleStiff(new double[] {E, Nu}, BVPType.PLANESTRAIN, K);
	    fs2.assembleMass(new double[] {density}, BVPType.CONSIST_MASS, M);
	    ScalarFunc myFunc = (coord, llabel, param) -> coord[1] > 0.99999? 1.0:0.0;
	    double[] x0 = fs.valueFromFunc(myFunc, null);
	    ScalarFunc betaFunc = (xy, llabel, param) -> beta;
	    
	    double[][] vecX = new double[1][];
	    vecX[0] = x0;
	    fs2.assembleHeatLoad(new ScalarFunc[] {betaFunc}, null, vecX, BVPType.COMMON, heatLoad);
	    
	    //����õ���ʼ���ٶ�
	    NewIterSSolver solverM = new NewIterSSolver(M);
	    solverM.PCGSSOR(heatLoad, acc, 1.5, 1);
	   
	    for(int i = 0, n = fs2.GetNdof(); i < n; i++) {
	    	acc0[i] = acc[i];
	    }

	    double[][] uv = new double[2][fs.GetNdof()];
	    fs2.extractComponent(acc0, uv);
	    mesh2D.toTecplot(dir + "acc.dat", uv);
	    MVO.fill(uv, 0.0);
	    
	    //K.sort();
	    //System.out.println(K);
	    //System.exit(0);
	    
	    SRMatrix ELMatrix = K.add((1.0 / (0.25 * dt * dt)), M);
	    //ELMatrix.sort();
	    //ELMatrix.toString();
	    //System.out.println(ELMatrix);
	    //System.exit(0);
	    
	    
	    String heatoutputFile = dir + "Heat100_2.dat";
	    mesh2D.toTecplot(heatoutputFile, time, x);

	    String elasticOutputFile =  dir + "Deform100_2.dat";
	    mesh2D.toTecplotDeform(elasticOutputFile, 0.0, uv, 0.0);
	    //call mesh2D.exportMeshDeformTime(elasticOutputFile, time, x, x, 0.0D0)
	    
	    double[] lineX = new double[lineNum + 1];
	    double[] lineY = new double[lineNum + 1];
	    
	    Arrays.setAll(lineX, i -> i * dis);

	    String XYLineFile = dir + "Line100_2.dat";
	    Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);

	    //x = 1.0D0
	    String expansion = dir + "Expansion100_2.dat";
	    FileWriter fw = null;
	    try {
			fw = new FileWriter(expansion);
			fw.write(0.0 + "\t" + 0.0 + "\n");
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	    //open(unit=20, file=expansion) 

	    //temp00File = trim(dir)//'Temp000_100_2.dat'
	    //open(unit=40, file=temp00File)

	    //XYStressYFile = trim(dir)//'StressYLine100_2.dat'
	    //call exportXYTime(XYStressYFile, lineX, time, lineY)
	    
	    fs.applyBC_MBN(LMatrix, Direct.All, outerLabel);
	    
	    
	    fs2.applyBC_MBN(K, Direct.Y, bottomLabel);
		fs2.applyBC_MBN(K, Direct.X, leftLabel);
	    
	    fs2.applyBC_MBN(ELMatrix, Direct.Y, bottomLabel);
		fs2.applyBC_MBN(ELMatrix, Direct.X, leftLabel);
	    
	    //ELMatrix.sort();
	    //System.out.println(ELMatrix);
	    //System.exit(0);
	    //call fs.applyBC(RMatrix, RHS, heatLabel, 1.0D0, 1)

	    //call fs2.applyBCZero(K, dispLabel1, directX)
	    //call fs2.applyBCZero(K, dispLabel2, directY)

	    //vel = 0.0D0
	    
	    NewIterSSolver SolverLM = new NewIterSSolver(LMatrix);
	    //NewIterSSolver solverL = new NewIterSSolver(LMatrix);
	    double[] Md = new double[fs2.GetNdof()];
	    double[] MMd = Md.clone();
	    
	    NewIterSSolver solverEL = new NewIterSSolver(ELMatrix);
	    time = 0.0;
	    double coef = 1.0 / ((1 + Nu) * alpha);
	    while(time < endTime) {
	        time += dt;
	        System.out.println("time: " + time);
	        //label(1) = 3
	        
	        //RHS = RMatrix * x
	        //call fs.applyBC(LMatrix, RHS, label, 1.0D0, 1)
	        fs.applyBC_MBN(RHS, Direct.All, 1.0, outerLabel);
	        SolverLM.PCGSSOR(RHS, x, 1.5, 1);
	        mesh2D.toTecplot(heatoutputFile, time, x);
	        RMatrix.mul(x, RHS);
	        
//	        for(int i = 0; i < lineNum + 1; i++) {
//	        	lineY[i] = fs.getValue(x, new double[] {0.25, lineX[i]});
//	        }
//	        Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);
	        
	        //call fs.valueXY(x, 0.25D0, 0.0D0, val)
	        //write(40, "(2F30.10)") time, val
	        //call exportXYTime(XYLineFile, lineX, time, lineY)

	        Arrays.fill(heatLoad, 0.0);
	        vecX[0] = x;
	        fs2.assembleHeatLoad(new ScalarFunc[] {betaFunc}, null, vecX, BVPType.COMMON, heatLoad);
	        //call fs2.assembleHeatLoad(beta, x, 0, HeatLoad)
	        //M.mul(1.0 / (0.25 * dt * dt), d, Md);
	        for(int ii = 0, nn = fs2.GetNdof(); ii < nn; ii++) {
	        	Md[ii] = 1.0 / (0.25 * dt * dt) * d[ii] + (1.0 / (0.25 * dt)) * vel[ii] + (0.5 / 0.25 - 1.0) * acc[ii];
	        }
	        M.mul(Md, MMd);
	        for(int ii = 0, nn = fs2.GetNdof(); ii < nn; ii++) {
	        	totalLoad[ii] = heatLoad[ii] + MMd[ii];
	        }
	        //System.out.println(MVO.toString(totalLoad));
	        //System.out.println(MVO.max(totalLoad) + "\t" + MVO.min(totalLoad));
	        
		    //fs2.applyBCZero(totalLoad, Direct.X, dispLabel1);
		    //fs2.applyBCZero(totalLoad, Direct.Y, dispLabel2);
	        
	        fs2.applyBC_MBN(totalLoad, Direct.Y, 0.0, bottomLabel);
			fs2.applyBC_MBN(totalLoad, Direct.X, 0.0, leftLabel);
			
		    solverEL.PCGSSOR(totalLoad, d, 1.5, 1);
		    fs2.extractComponent(d, uv);
		    mesh2D.toTecplotDeform(elasticOutputFile, time, uv, coef);
		    
	        //call fs2.computeStrain(u, v, ex, ey, exy)
	        //call fs2.computeStress(ex, ey, exy, u, E, Nu, tx, ty, txy, u, 1)

	        //do i = 1, lineNum + 1
	        //    call fs.valueXY(ty, 0.25D0, lineX(i), val)
	        //    lineY(i) = val 
	        //end do
	        //call exportXYTime(XYStressYFile, lineX, time, lineY)

		    for(int j = 0; j < fs2.GetNdof(); j++) {
		    	acc[j] = 1/(omega*dt*dt) * (d[j] - d0[j]) -1/(omega * dt) * vel0[j] - (0.5 /omega - 1) * acc0[j];
		    	vel[j] = vel0[j] + ((1-0.5) * acc0[j] + 0.5 * acc[j]) * dt;
		    	acc0[j] = acc[j];
		    	vel0[j] = vel[j];
		    	d0[j] = d[j];
		    }
	        
		    try {
				fw.write(time + "\t" + fs.getValue(uv[1], new double[] {0.0, 1.0}) * coef + "\n");
			} catch (Exception e) {
				e.printStackTrace();
			}
	        //call fs.valueXY(v, 0.25D0, 1.0D0, val)
	        //write(20, "(2F20.10)") time, val * (1-Nu)/ (alpha * (1 + Nu))	    	
	    }
	    try {
	    	fw.flush();
			fw.close();
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	}

}
